Hydrodynamic modes in a confined granular fluid 
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Confined granular fluids, placed in a shallow box that is vibrated vertically, can achieve homo- 
geneous stationary states thanks to energy injection mechanisms that take place throughout the 
system. These states can be stable even at high densities and inelasticities allowing for a detailed 
analysis of the hydrodynamic modes that govern the dynamics of granular fluids. Analyzing the 
decay of the time correlation functions it is shown that there is a crossover between a quasielastic 
regime in which energy evolves as a slow mode, to a inelastic regime, with energy slaved to the other 
conserved fields. The two regimes have well differentiated transport properties and, in the inelastic 
regime, the dynamics can be described by a reduced hydrodynamics with modified longitudinal 
viscosity and sound speed. The crossover between the two regimes takes place at a wavevector that 
is proportional to the inelasticity. A two dimensional granular model, with collisions that mimic the 
energy transfers that take place in a confined system is studied by means of microscopic simulations. 
The results show excellent agreement with the theoretical framework and allows the validation of 
hydrodynamic-like models. 



o 

CO 



-a 
a 
o 

(N 
> 
00 



o 



x 



I. INTRODUCTION 

Granular fluids have become a prototype of non- 
equilibrium matter. The need of permanent energy in- 
jection to counter-balance the energy dissipation in the 
grain interactions, place these systems under permanent 
non-equilibrium conditions. Energy is injected through 
boundaries or external fields and it is dissipated at the 
small scale of grain-grain collisions. This fact violates 
the detailed balance condition necessary to reach equi- 
librium. It is then, one of the objectives in the study of 
granular fluids, the construction of valid statistical me- 
chanics tools under these non-equilibrium conditions 

The usual approaches to the study of granular fluids 
are kinetic theory (with different levels of approximation) 
or hydrodynamic-like models for the relevant fields Q. 
The dissipative nature of collisions implies that granular 
media cannot be simultaneously in homogeneous and sta- 
tionary states and, typically, spatio-temporal structures 
develop [1, The homogeneous cooling state (HCS), 
in which the energy is non-stationary, has been widely 
studied showing that it becomes unstable in the long- 
wavelength regime It is the reference state for 
developing kinetic and hydrodynamic models of granu- 
lar media with small or vanishing driving [Tol - [T2j . Con- 
versely, stationary states can only be obtained by perma- 
nent energy influx. When granular media are driven by 
boundaries, typically large inhomogeneities develop even 
in the stationary regimes (see for example, [l3[ ). Local 
energy balance can be obtained by compensating the en- 
ergy dissipation with shear heating. In this case it is of 
particular interest the uniform shear flow (USF), in which 
all fields are uniform except for the velocity that shows a 
linear profile [3 EH- As in the HCS case, the USF serves 
as a reference state to develop kinetic and hydrodynamic 
models. An important outcome is that the transport co- 
efficients for the linear dynamics close to the HCS and 



USF states are different [3, [l5[. In both cases and in 
other studied states, however, a generic feature appears. 
The evolution of the energy shows two well differentiated 
regimes, depending on the dissipation (T^]. At low dissi- 
pations the energy evolves in long time scales and can be 
treated as another hydrodynamic field in equal foot as 
the conserved fields (density and momentum). At large 
dissipations, on the other hand, the energy evolves fast 
and it is slaved to the density and velocity field. An 
example of this slaving is found in avalanches, in which 
the granular temperature is proportional to the vel ocity 
gradient squared in the so-called Bagnold scaling (l7lll8|. 

The crossover between the previous regimes, the 
quasielastic and the inelastic ones, is difficult to observe 
and characterize qualitatively. Only under dilute condi- 
tions, the quasielastic regime is observable at finite in- 
elasticities. At moderate densities the inelasticity must 
be extremely small otherwise the only visible regime is 
the inelastic one [lg, [l9( . Related to this is the fact that 
in dense or moderately dense regimes, granular fluids de- 
velop large inhomogeneities, when the use of hydrody- 
namic equations (with transport laws linear in the field 
gradients) are of questionable validity 0. Some other 
authors extended the hydrodynamic description by us- 
ing nonlinear constitutive relations [2(J HH • Homogene- 
ity can be achieved in small systems, in which the un- 
stable wavevectors are not accessible . Again, the lim- 
itation to large wave vectors renders hydrodynamics of 
limit validity. In summary, the regime crossover has not 
been tested under dense inelastic conditions, issue that 
is studied here. 

The quasi two dimensional (Q2D) geometry offers a 
possibility to study this crossover and the properties of 
the hydrodynamic modes near stationary and homoge- 
neous regimes. In this geometry, grains are placed in a 
box with large horizontal dimensions, while the vertical 
dimension is small, typically less that two diameters in 
height. When the box is vertically vibrated, grains get 
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energy through the collisions with the top and bottom 
walls and this energy is then transferred to the horizon- 
tal degrees of freedom via grain-grain collisions. As these 
collisions are also inelastic, the system can achieve sta- 
tionary states with finite energies. The vertical scale is 
fast and evolves in the scale of a few vibration periods. 
The horizontal dynamics, on the other hand, evolves in 
larger times scales characterized by the density and mo- 
mentum conservation. In this geometry, it is known that 
in a wide range of parameters including dense inelastic 
conditions, the system remains homogeneous in the hor- 
izontal directions (22| - [25j . The key element that allows 
for the establishment of stationary homogeneous states is 
that, for the effective horizontal dynamics, there is a dis- 
tributed energy injection source. In the case of the Q2D 
systems, in the absence of friction, this energy source is 
Galilean invariant and conserves momentum locally. 

In this article we study the hydrodynamic modes in a 
granular fluid with a distributed energy injection mecha- 
nism similar to the one in the Q2D geometry. The anal- 
ysis, although inspired in the Q2D geometry, is generic 
and valid for three dimensional systems if a distributed 
energy injection mechanism is devised. It will be shown 
that there is a crossover between the quasielastic and in- 
elastic regimes and the properties of the modes will be 
studied in detail in both regimes. The analysis will be 
done studying the density-density correlation functions 
that are obtained from fluctuating hydrodynamics. The 
intermediate scattering function and the dynamic struc- 
ture factors provide information of the relevant modes 
and their time dependence. Finally, we present a dis- 
crete microscopic model in which grains can gain or dis- 
sipate energy at collisions and the results obtained from 
molecular dynamics simulation are analyzed under the 
described framework. 

The plan of the paper is the following. In Sect. [TT| we 
develop the framework for the analysis of the hydrody- 
namic modes using correlation functions for a granular 
fluid. In Sections Mil and HVl the inelastic and quasielas- 
tic regimes and analysed in detail. In the inelastic case, 
we derive the temperature slaving that gives rise to a re- 
duced hydrodynamics. Section IVl analyzes the crossover 
between these regimes, showing that it takes place in 
a wavevector proportional to the inelasticity. A micro- 
scopic collisional model in two dimensions, that mimics 
the Q2D dynamics is presented in Sect. IVII Simulations 
of this model and comparison with the theoretical frame- 
work are shown in Sect. IVIII Finally, conclusions are 
given in Sect. IVIIII 



II. HYDRODYNAMIC MODES FOR 
GRANULAR FLUIDS 

The goal of this section is to follow a procedure equiv- 
alent to that of Landau and Placzeck for granular flu- 
ids, considering the particular issues of such systems, like 
modification of the hydrodynamic equations, and deriv- 



ing the modifications in both the intermediate scatter- 
ing function F(k,t) and the dynamic structure factors 
S(k,u>). Time correlation functions of equilibrium fluc- 
tuations are standard tools in the study on fluids, as they 
contain equilibrium properties (like, e.g. specific heats, 
or the speed of sound) as well as non equilibrium ones 
(transport coefficients). Onsager's regression hypothesis 
states that spontaneous fluctuations in equilibrium obey 
the same evolution equations that describe the macro- 
scopic relaxation of an external perturbation, provided 
that the perturbation is weak. As, in the hydrodynamic 
limit of long wave lengths, macroscopic relaxation pro- 
ceeds via the Navier Stokes equations, the correlation 
function also evolves according to those equations. Then, 
the correlation functions can be used to measure trans- 
port coefficients and other thermodynamical properties 
of the fluid [H,[13. 

Let us define a general space and time correlation func- 
tion between the dynamic variables A and B as 

C AB {v,t) = {5A{v + v',t + t')5B{v'X)), (1) 

where 5A(r,t) = A(r,t) — (A(r,t)} is the fluctuation of 
the variable A with respect to its average value. In the 
definition above, we have assumed that the system is spa- 
tially homogeneous and invariant under time translation, 
so the system must reach a stationary state; otherwise the 
correlation function would depend on r and r' and also 
on t and t' . Although the definition for Cab is general, 
we will restrict to the study of density autocorrelation 
function, where A = B = p, or the velocity correlation 
function, where A = B = u, that leads to a tensorial cor- 
relation function. Definitions of the observables in term 
of microscopic quantities are given in Appendix |A"1 

For practical purposes, it is convenient to take the 
Fourier transform in space of CAs(r,£) to obtain the 
so called intermediate scattering function, denoted by 
F(k,t). Furthermore, the Fourier transform in time can 
be taken to get the dynamic structure factor, S(k,uj), 
whose properties for equilibrium fluids are explained 
in Appendix [A] At long or hydrodynamic wavelengths, 
larger than the mean free path or the size of the 
molecules, the so called Landau-Plazceck approximation 
allows calculation of F(k,t) and S(k,uj). Such calcula- 
tion is based in the fact that the time dependent hydrody- 
namic fields can be described by the set of Navier Stokes 
equations at linear level. In contrast with the molecular 
fluids, S(k,ui) has not been widely used for inelastic or 
dissipative systems, but only recently [l9L [28L [29l| . 

In order to construct such functions, we need the evo- 
lution equations for the system. We consider granular 
particles with energy is dissipated at every collision, so 
we will not consider here systems with a Stokes friction, 
like those of [28, 29]. Then, in order to reach a stationary 
state, we have to supply energy into the system. There 
are many models for energy injection [ijj [28| - |3^ |. and 
we will introduce a collisional model for energy injection 
in Sect. IVII but for the time being we will develop the 
theory as much as we can without specifying its detailed 
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form. We will assume that the thermostat does not inject 
momentum, but only energy, and it is Galilean invariant. 
Under such assumption, the equations for the density 
field and momentum density are those of usual fluids: 
the continuity equation and the Navier Stokes' one. 

Despite the equation for the energy (or granular tem- 
perature T) does not derive from a microscopic conserved 
quantity, one can write a balance equation for it. Taking 
as a starting point the conservation equation for the en- 
ergy for elastic fluids, we must add a term that accounts 
for the dissipation and the energy injection. We will de- 
note such term in the temperature equation by G(p,T), 
where we make explicit the dependence on the density 
and the temperature. It will also depend on microscopic 
coefficients like, for instance, the coefficient of normal 
restitution a, and also on parameters that characterize 
the energy injection. 

The second modification comes from the constitutive 
relation for the energy flux [3j, llfj, |33j, |34| . It includes the 
usual heat conduction term, given by Fourier's law and 
besides, there is a contribution proportional to the gradi- 
ent of the density, that has no counterpart in molecular 
fluids. Then the heat flux reads 



-kVT — /iVp. 



(2) 



With such considerations into account, we can write 
the nonlinear temperature equation as 

d t T(r,t) = -u.VT- — (Q) V u 
pc v \dTJ p 

— — : Vu+ — V(kVT + /iV/3) - G(p,T). (3) 
pc v pc v 

Here cv is the specific heat at constant volume, p is the 
hydrostatic pressure and P' is the traceless part of the 
stress tensor. At the hydrodynamic level, we assume that 
the stress tensor is Newtonian, characterized by shear 
and bulk viscosities rj and r\v ■ The arbitrary minus sign 
in front of G has been included for later convenience. 

Let us note again that this equation does not derive 
from a microscopical conserved quantity, reflecting that 
the term that describes the dissipation and the energy 
injection, G, does not derive from a flux term, and there- 
fore is not proportional to a gradient. 

When the system of granular particles evolves, the 
temperature may reach a stationary value, which is a 
balance between the dissipation and the energy injection. 
We assume that there exists a homogeneous stationary 
state. We can calculate the stationary temperature T st 
by integrating over the whole system Eq. ([3]), where all 
terms under a spatial derivative vanish, arriving at the 
expression 



G(p,T st ) =0. 



(4) 



This equation defines the stationary temperature in 
terms of the density, and other parameters included in 
G, like dissipation or the energy injection, that sets the 
functional form of G. 



When studying fluctuations about the stationary state, 
we use Onsager's regression hypothesis. We linearize the 
evolution equations around the stationary state charac- 
terized by a constant density, a vanishing velocity and the 
temperature T st . Then, we define fluctuations around 
such state as 



p(r,t) = p + Sp(r,t), 
u(r,t) = Su(r,t), 
T(r,t) =T st + ST(r,t). 



(5) 
(6) 
(7) 



Linearization around such steady state follows the usual 
procedure as for molecular fluids. The new term, G(p, T) 
linearizes to first order as 



G(p, T) ~ G p Sp + G T ST, 



(8) 



where Gx denotes the derivative of G respect to the vari- 
able X evaluated at the average density and the station- 
ary temperature. In an elastic fluid, the terms Gx are 
absent. They are present only in dissipative media and, 
in fact, they are proportional to the inelasticity of the 
medium. In what follows, before giving any explicit form 
of the energy injection mechanism, we will refer to them 
as the dissipation terms. 

Then, the set of linear equations in the Fourier vari- 
able V — > ik reads dt^f = — M^P. The vector W contains 
the Fourier transform of the fields, and M is the pseudo 
hydrodynamic matrix, with expressions 



6p(k,t) 
5u\\(k,t) 
Su±(k, t) 

ST{k,t) 



ikp p 
P 





G 



k z /i 



ikp 


ikT at pq 





k 2 v 






ikpT 
P 





G T + — 



(9) 



(10) 



where px denotes the derivative of p respect to X eval- 
uated at the average density and the stationary tem- 
perature, v = r\j p is the kinematic viscosity and vi = 
(jl + Vv)/p is the longitudinal kinematic viscosity. As 
usual, we have decomposed the velocity field u(k, t) into 
its longitudinal, uii(k, t) = k- u(k, t), and transversal, 

u±(k,t) = u(k, t) — ktiy(k,t), parts. The transversal 
part is in fact, a D — 1 dimensional vector, and so it is 
the matrix that contains their components. 

The matrix M is the modified hydrodynamic matrix 
for granular fluids. It differs from the hydrodynamic ma- 
trix for molecular fluids in two elements, related with the 
temperature: the (T, p)-element includes the new trans- 
port coefficient p coming from Eq. ([2]l. and the term G p , 
while the (T, T)-element contains the term Gt- Such 
terms modify drastically the spectrum of the matrix M. 
The matrix also has the modified equation of state and 
transport coefficients for a granular fluid, but these only 
modify quantitatively the matrix elements. 
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Obtention of the time dependence of the fields, re- 
quired in Eq. (JlJ , involves the diagonalization of the ma- 
trix M . As the matrix is not Hermitian there are two sets 
of orthonormal eigenvectors, right and left ones, given by 

Mipi = Xiipr, (f>M — Xi4>i\ 4>i ■ ipj — Sij, (11) 

with components labelled by the superindex /3, that can 
take the values: f3 = (p, ||,_L,T), in a self-explanatory 
notation. Then, the solution for the deviations at time i, 
denoted by ^(i) are 

*(t) = ^e- Ai ViCi, (12) 
i 

where the coefficients are the projections of the fluc- 
tuations at initial time over the left eigenvectors 

c i =<j) i - V(t = 0) 

= tfSpQt) + 4S|(k) + <j>t5u x (k) + <f>f5T(k), (13) 

and the fluctuations without explicit dependence on time 
are evaluated at t — 0. As we see, all the time dependence 
is contained in the exponential terms, of those the Fourier 
transform will be taken to get S(k,uj). 

As we are mainly interested in the density-density cor- 
relation function (see, however, Sec. IVIFDl when we also 
study the transversal correlation function), we need to 
multiply the density fluctuation at time t with that at 
time zero, obtaining 

F(k,t) = ^-{6p(k,t)6p(-k,0)) (14) 
pV 

where V is the volume. Here S/3 p (k) is the static (equal 
time) structure factors between the held /3 and density 
held p. In equilibrium, such structure factors are diag- 
onal, that is, only the density-density term, S pp , con- 
tributes to the sum in Eq. (fT4"]l [271] . However, in non- 
equilibrium fluids, the structure factors are not diagonal 
[35j . For symmetry reasons scalar and vectorial fields do 
not couple, implying that S\\ p = S± p = 0, while S p t ^= 0. 



Such static structure factors can be calculated for gran- 
ular fluids by using, e.g. the technique developed in [3(| 
for a 'random kick' driving. 

Before doing the full diagonalization, we note that M 
is positive definite for small wavevectors, that is all hy- 
drodynamic modes are stable, as long as GtP p > G p pt, 
otherwise one mode becomes unstable. This instability is 
of van der Waals type, related to the negative compress- 
ibility of the reduced dynamics at small wavevectors (see 
the end of Sec. HITO [H HI- In Sec. ED we will show 
that this condition is always fulfilled for the energy in- 
jection method we devise. 

From the structure of the pseudo-hydrodynamic ma- 
trix, the transverse mode decouples from the rest obtain- 
ing directly the associate eigenvalue. The other three 
modes couple and give contributions to the dynamic 
structure factor. Considering the parity and complex 
structure of M it is possible to deduce that the eigen- 
values have the form 

A± = ±iuj B {k) + f(fc), (16) 
A T = D T (k), (17) 
A_l = ffc 2 , (18) 

where A± are the eigenvalues associated with the sound 
modes and At to the heat mode. Using the standard 
notation of elastic fluids Dt (even function in k) is the 
dissipation rate of the thermal modes, V (even in k) is 
the dissipation rate of the sound mode, and ujb (odd in 
k) is the frequency of the sound modes. 

In the Landau-Plazceck theory of elastic fluids, Dt = 
Dxk 2 , r = Tk 2 and uib = c s k, where Dt is the thermal 
diffusivity, T is the sound damping constant, and c s is the 
adiabatic sound velocity. In next sections, the eigenval- 
ues of the inelastic model are computed and two regimes 
are clearly differentiated: the so called dissipative regime 
(Sect. IIIip . usual in granular media, and the quasielastic 
regime (Sect. IIV|) . The latter deals with the elastic limit, 
to make connection with the usual hydrodynamics, and 
to verify that the elastic limit is a singular limit. The 
crossover is studied in Sect. [V] Finally, comparison with 
molecular dynamics simulations of a collisional model is 
done. 

Starting from Eq. (I14[) and considering the temporal 
parity of F(k, t) and the presence of the three hydrody- 
namic modes that couple to the density, the intermediate 
scattering function can be written as 



F(k,t) = S(k) 



(1 - 7- 1 ) e- B - 1 + ( 7' 1 cos( Ws t) + F+( 7 , 1)Dt sinCwflt) ) e~ ft 



(19) 
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and, therefore, the dynamic structure factor is 



S(k,u) =S(k) 



2(l- 7 -) 



D 2 T 



uj 2 + D\ 



+ 1 



r 



r + ( 7 -i) J D T 



(lo + lj b ) 2 +T 2 {lo-u b ) 2 + T 2 , 

CO — LOB 



l^B 



(lo + lubY+T 2 {lj-lub) 2 + T 2 



where S(k) is the static (density-density) structure fac- 
tor. These expressions are formally equal to their coun- 
terparts in elastic fluids (26|, containing the Rayleigh, 
Brillouin and asymmetric peaks. However, the factor 7, 
that in equilibrium is the adiabatic constant or the ratio 
of specific heats, contains here the static structure fac- 
tor St p - Geometrically, 7 is defined such that the area 
enclosed by the thermal peak (Rayleigh) divided by the 
area under the sound peaks (Brillouin) is 7 — 1 . 



A. Energy scaling and dimensionless variables 

To simplify the analysis, we consider the case of in- 
elastic hard particles of diameter a and mass m, charac- 
terized by a velocity-independent restitution coefficient. 
Before giving details of the specific implementation of the 
energy injection mechanism, we assume that it introduces 
a unique energy scale. Fields and time can be rescaled 
according to this energy scale. Instead we use the equiv- 
alent procedure of rescaling to the stationary tempera- 
ture T st . Therefore, the fields, transport coefficients, and 
eigenvalues rescale by dimensional arguments as: p — > p, 
u -> VT^u, T -> T st f, v -> VT^V, vi -> VWv u 
K -> VT^K,J^ -> (T st f/ 2 p, p T ->• p?, Pp -> T s %, 
G T -> VT^Gt, G p -> {T st ) 3 / 2 Gp, and A -> y/T^X. 
Finally, dimensionless magnitudes are defined by further 
normalization with appropriate powers of m and a for 
each quantity. To simplify notation, the hat symbol will 
be suppressed in what follows. 



III. INELASTIC REGIME 

The first regime we consider is the inelastic regime, 
in which dissipation is large enough such that the en- 
ergy evolves faster than the density and momentum at 
the relevant, small, wavevectors. The eigenvalues in this 
regime are obtained as a series in k keeping finite the 
dissipation parameters G p and Gt- The results, up to k 2 
(i.e. hydrodynamic) order, are 



*-* + (£-%r-^)* (21) 

A ± = vk 2 . (22) 

The thermal mode is drastically modified compared to 
elastic fluids, as contains a zero order contribution be- 
cause of lack of energy conservation. Therefore, it is not 
a slow mode anymore. 

The sound mode contains the dissipative sound speed, 
c c i = y/pp — ptG p /Gt, instead of the adiabatic veloc- 
ity like in elastic systems (see Appendix |BJ . Again, it is 
a consequence of the lack of energy conservation. The 
dissipative sound velocity corresponds to the isothermal 
velocity ^/Pp', corrected by the instantaneous coupling be- 
tween T and p given by the constraint (j4]). The sound 
modes are stable as long as p p Gt > PtG p as was pre- 
viously indicated in Section [TTJ The other modes are un- 
conditionally stable for small wavevectors. The viscous 
mode is not coupled to the dissipation terms and there- 
fore is the same as in the elastic case. 

Note that the sound velocity, the sound damping con- 
stant and the k 2 term in the thermal mode do not ap- 
proach the elastic values when the dissipation vanishes 
(G p — > and Gt 0). The elastic limit is singu- 
lar and a detailed scaling must be performed to match 
both regimes. This scaling is done in Section IIVI and 
the numerical demonstration of the crossover between the 
regimes is done in Section fVIIl 

Finally, let us mention that the new transport coeffi- 
cient p does not appear in the eigenvalues (|2T)|) - (|2"TT) of the 
matrix M, up to order k 2 . Therefore, one cannot measure 
such transport coefficient with the method developed in 
this paper, but has to devise different methods |34j |. 

A. Slaving of the temperature field 

Being the temperature a fast field, it is possible to slave 
it to the density and velocity fields and obtain a simpler 
dynamics. The linearized temperature equation reads 

/„ k 2 u\ „ ikpx „ /„ k 2 n\ „ 

d t ST = - [Gp + — Sp = —5m - [Gt + ) ST. 

\ c v pj c v p \ c v pj 

(23) 
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At the time scale of the slow modes the temperature 
has relaxed to the stationary solution of (f2"3"|) allow- 
ing to obtain an explicit expression for the temperature 
field. As ST enters as a spatial derivative in the mo- 
mentum equation ([9J only terms up to order k must 
be retained to reproduce hydrodynamic modes, giving 
$T = ~73^3p ~ cvpG T " However, when substituted in 
the momentum equation the resulting eigenvalues differ 
to the sound modes ([20]) in terms of order k 2 . The reason 
is that this simple Markovian slaving does not consider 
all contributions to the same order in k. Formally, equa- 
tion (l23l) can be written as 



d t ST= -jST+S(t), 



(24) 



with S = — (Gp - 
{G T + k 2 n/{c vP )y 



ST(t) 



k 2 p/(c v p)) 5p - % -^jSu\\ and 7 = 
. Its solution is 



dt' e-^'S{t-t') 



dt' e~ 7 *' 



S(t) S(t) 



S(t) - S(t)t' 



7 



r 



(25) 



where the first term gives the Markovian slaving. To 
evaluate the second term, the hydrodynamic equations 
for Sp and 8u\\ are used, noting that only terms up to 
order k must be retained. A third term, proportional to 
S gives contributions of order k 2 , not being relevant to 
the hydrodynamic model. The resulting slaving for ST is 



ST 



G p ^ ikp T ^ 
Gt cypGr " 



ikpGp 
G 2 ^ 



Sun 



(26) 



and the reduced hydrodynamic matrix for the slow vari- 
' S P (k,t) 
<Wk,i) 



ables * 



M = 







ikp 



\ p G T p J \ Gfp cVGtP' 1 J 



(27) 



IV. QUASIELASTIC LIMIT 

When dissipation is small it is expected that the tem- 
perature field becomes a slow field together with the den- 
sity and momentum fields. The dissipation terms Gx 
must be compared with the terms proportional to k 2 in 
(f2"Tj) . Therefore, the scaling that capture this quasielastic 
regime must be done in the dissipation together with the 



wavevector. It is obtained doing G p 



e Gp, Gt 



e 2 G 7 



and k = ek (e is a formal small parameter) and computing 
the eigenvalues as series in e. Keeping terms up to order 
e 2 and going back to the original variables the eigenvalues 
are 

A + = PT{G TPT /cyp 2 + Gp) ± ^ + ^ 



2c 2 . 



At = 



(G T p P - G p p 2 T ) 



D T k 2 , 



A 1 =vk 2 



where 



\lpp+p T /c v p 2 



r 

D T 



2 2c v p{p T +p p p 2 c v ) 

KPpP 



p T + c v ppp 2 



(29) 
(30) 

(31) 
(32) 
(33) 



are the elastic adiabatic sound speed, sound damping 
constant and the thermal diffusivity, respectively. 

The fc-independent term in the sound modes could lead 
to the erroneous impression that they are fast modes and 
do not correspond to conserved fields. This is not the 
case as it should be recalled that the quasielastic regime 
is obtained when the dissipation scales as k 2 and there- 
fore the fc-independent terms vanish also when the wave 
vector goes to zero. A detailed analysis of the crossover 
of the two regimes and a geometric interpretation of the 
scaling is given in the next section. 

Note that, as in the inelastic regime, the eigenvalues 
do not depend on the transport coefficient /x, up to order 
k 2 . 



Using this reduced hydrodynamic model, the sound 
modes are correctly reobtained. Note that in the reduced 
hydrodynamics, there is an effective longitudinal viscos- 
ity that is the bare longitudinal viscosity modified by 
a term that depends on the energy injection-dissipation 
mechanism. 

In this reduced dynamics, it is clear the effect of the 
stability condition p p Gt > PtG p . It guarantees that 
the sound modes are stable. When it is not fulfilled, 
a spinodal decomposition develops and non-linear terms 
are necessary to describe the long term dynamics, result- 
ing in the van der Waals normal form (3q . l37j . 



V. CROSSOVER WAVEVECTOR 

The two regimes described in the previous Sections (in- 
elastic and quasielastic) lead to different eigenvalues of 
the hydrodynamic modes. We make a special remark on 
the sound velocity which is either isothermal or adiabatic. 
For a given dissipation, there is a crossover wavevector 
ko, such that when k -C fco the dynamics is inelastic while 
if k 3> fco the dynamics is quasielastic (see Fig. [T]) . 

To compute fc , we consider the full expression of the 
eigenvalues associated with the sound mode and get the 
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real part of the eigenvalues, T. For a given small dissipa- 
tion, kg is the wavevector in which T has changed appre- 
ciably from the k — limit. It can be verified by doing 
a full diagonilization of the pseudo-hydrodynamic matrix 
that, in the limit of small dissipations, T ps Tq— T4fc 4 +. . .. 
Therefore, the crossover wavevector is computed as 



k 4 = - 

K — 



24r\ 

9 4 f J 



(34) 



In the limit of small dissipation (i.e. the dimensionless 
non-hydrodynamic terms Gt and G p are small) ko oc Gx ■ 
The crossover wavevector is proportional to the dissipa- 
tion with a proportionality constant that depends on the 
equation of state, the heat capacity and the ratio G p /Gt- 
The regimes are presented schematically in Fig. [T] In 
the elastic limit, where Gx — > 0, only the quasielastic 
regime, in the appropriate limit, is a valid description. 
When the dissipation is finite, the relevant regime de- 
pends on k, obtaining the inelastic regime in the limit of 
large systems (k going to zero) . The figure also shows the 
singular character of the elastic limit: it is not possible 
to obtain it for finite dissipations by taking k —± 0. It 
can only be obtained if the inelasticity is reduced simul- 
taneously as shown in the scaling of Section IIVI 



G X k 



Inelastic 
regime 



Quasielastic 
regime 



k 



FIG. 1. Schematic representation of the quasielastic and in- 
elastic regimes, in terms of the wavevector k and dissipation 
Gx- The crossover wavelength is ko oc Gx- At finite dissipa- 
tion, both regimes are present, depending on the wavevector. 
Only in the elastic case, the quasielastic regime is valid for all 
wavevectors. 



VI. COLLISIONAL MODEL 

A particular geometry that has gained interest in the 
study of granular media is the quasi two-dimensional one 
(Q2D). Here, the box is large in the horizontal directions, 
while the vertical one is smaller than two particles' di- 
ameter, such that grains cannot be on top of another. If 
the box is vertically vibrated, with a maximum acceler- 
ation larger than gravity, grains gain vertical energy by 



collisions with the top and bottom walls and this energy 
is transferred to the horizontal directions through grain- 
grain collisions. Seeing from above the granular system 
is fluidized and can remain homogeneous under a large 
range of parameters. Varying the vibration amplitude 
and frequency the system develops a phase transition 
mediated by waves |25l| with a solid-like region coexisting 
with the fluid 22, 38] . Here we focus on the homogeneous 
state and with that purpose an effective two-dimensional 
model is proposed. 

If only the horizontal two-dimensional degrees of free- 
dom are considered, collisions can either dissipate or gain 
energy, depending on the amount of vertical energy grains 
have and the restitution coefficients. This idea was ex- 
ploited in Ref. (32j in which the restitution coefficient 
was a random variable with possible outcomes larger than 
one. That model, however, lacked of an energy scale and 
the total energy of the system performs a random walk, 
not reaching a steady state. In the Q2D system, the ver- 
tical energy scale of the grains is fixed by the vibration 
parameters and so is the typical energy that is transferred 
from the vertical to the horizontal degrees of freedom. 
We propose a two-dimensional hard disk model, in which 
collisions are characterized by a constant restitution co- 
efficient a and an extra velocity A that is is added to 
the relative motion. This extra velocity points outwards 
in the normal direction a as required by conservation of 
angular momentum [39] . The collision rule for the post- 
collisional velocities reads 



1 



v* = vi - -(1 + a)(vi2 • o)o - a A, 



(35) 



v 2 = v 2 



(I + a)(v 12 ■ o)a + a A, 



where v 12 = Vj — v 2 is the relative velocity, a points from 
particle 1 to 2, and particles are approaching if vi 2 -<r > 0. 

With this set of collision rules, momentum is con- 
served, but energy is not conserved. The energy change 
in a given collision is 



E* -E 



mA 2 + m(vi 2 • a)aA — m(vi2 ■ cr) 



(36) 

I -a? 

4 ' 
(37) 



Considering a Maxwellian velocity distribution, absence 
of velocity correlations and static pair correlation func- 
tion at contact x> the energy dissipation rate per unit 
area, that should be included in the hydrodynamic equa- 
tions, is 



G = - 



L0( P ,T) 



iA 2 + aAvWT - T(l - a 2 ) , (38) 



where u(p,T) — 2pax^/^T/m is the collision frequency 
and the prefactor 1/2 compensates the double counting 
of collisions. 

The resulting expression of G has the remarkable prop- 
erty that it is factorized into two terms that depend only 
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on p and T, respectively. This feature is a result of en- 
ergy being injected and dissipated at collisions but not 
on the particular way of the collision rule (|35|) . As a 
consequence, the stationary temperature, T st , is density 
independent and is given by 

The stationary temperature diverges in the elastic limit 
(a — > 1) as energy is injected in every collision but no 
dissipation takes place. This divergency is correctly re- 
produced in simulations (see Sect. IVTIj) . In order to ob- 
tain a finite temperature in the elastic limit, both the 
dissipation and the energy injection must vanish (a — > 1, 
A -> 0, and A/(l - a) -> const.). 

In hard sphere models, the pressure is the tempera- 
ture times a monotonic function of the density. As in the 
proposed model (thermostated at collisions) the station- 
ary temperature is density independent, the pressure in- 
creases monotonically with density and no effective neg- 
ative compressibility can be produced (37j . Therefore, 
effective two-dimensional collisional models are not able 
to reproduce the solid-liquid transition in Q2D systems. 
Our event driven simulations show that indeed the sys- 
tem in stable for all densities and inelasticities, giving 
rise to stationary homogeneous states. 

The factorization of G also implies that G p = and 
therefore many expressions of the previous sections sim- 
plify. In particular the dissipative sound speed reduces to 
c d = yfPpi t ne isothermal sound speed. Also, the hydro- 
dynamic matrix is positive definitive for any parameters 
and all modes are stable. 

Finally, as advanced in Sect. [11] the dissipation terms 
are proportional to the inelasticity. Indeed, when the par- 
tial derivative Gt is computed and the stationary tem- 
perature is substituted, its dimensionless form is propor- 
tional to the inelasticity 1 — a. 



VII. NUMERICAL RESULTS OF THE 
COLLISIONAL MODEL 

The effective 2D collisional model is simulated using 
the event driven algorithm for hard disks, considering 
the collision rule (|35|) . The disk diameter a, particle mass 
to, and the energy injection parameter A are used to fix 
length, mass and time units. Simulations are done for 
systems of different restitution coefficients a and global 
density p. In this case we compute directly the interme- 
diate scattering function, F(k,t), that for k ^ reads 

F ^ = ^X>PHk- (r,(t) - r,(0))^ . (40) 

For isotropic steady states, as those reached for our 
model, F(k,t) — F(k,t), where k = |k|. The reason 



to compute F(k, t) instead of C(r, t) is that such expres- 
sion of the intermediate structure factor is well suited for 
numerical calculations, as it deals with analytical func- 
tions (exponentials) of the positions of the particles. Such 
exponentials have to be calculated at regular instants of 
time t for several k vectors. Then, time correlations must 
be performed to obtain F(k,t). On the contrary, the 
real space correlation functions, defined in Eq. ([T|) either 
require handling a delta function, or either performing 
averages over certain spatial domains in order to obtain 
coarse grained densities. Once F(k,t) is calculated, we 
carry out numerical Fourier transforms in time to com- 
pute S(k,uj). To explore different wavevectors, different 
box sizes (and number of grains N) were used to increase 
the wavelength and also, different wavevectors were an- 
alyzed simultaneously for a given system size. Finally, 
different aspect ratio L x /L y were explored. In all cases, 
the system was verified to be stable and statistically ho- 
mogeneous and, for a given wavevector, the computed 
physical quantities are independent of system size or as- 
pect ratio. 

The systems are initialized with a homogeneous dis- 
tribution in space and velocities are sorted according to 
a Maxwellian distribution at the theoretical temperature 
T st (|39|) . Then, the system is let to relax until a sta- 
tionary state is reached. Figure [5] shows the station- 
ary temperature obtained in simulations compared with 
the predicted value. The agreement is excellent even for 
large inelasticities (a close to 0), where the hypothesis 
of Maxwellian distribution or absence of velocity corre- 
lations are expected to fail. Note that there is a small 
density dependence (few percents) on the stationary tem- 
perature implying that either there are velocity correla- 
tions or the velocity distribution depends on density. 

In the stationary state, the density fluctuations are 
obtained for different wavevectors, computing the in- 
termediate scattering function F(k,t) and the dynamic 
structure factor S(k,tu), that is shown in Fig. [3] for 
two different wavevectors. The two regimes (inelastic 
and quasielastic) are clearly seen: in the first case, only 
two peaks (sound modes, coming from Eq. (27)) are 
present while in the second case the three peaks (sound 
and heat modes) are visible. In all cases, no more 
than three peaks are observed showing that the pseudo- 
hydrodynamic model describes correctly the dynamics 
and it is not necessary to introduce additional kinetic 
modes [II [HI. 

In what follows we focus on the intermediate density 
p = Na 2 1 (L x L y ) — 0.4; similar results are obtained for 
other densities. To analyze the data, instead of fitting the 
position and width of the peaks in S(k,u), we directly 
fit F(k,t) with (TT5)) to obtain 7, lob, T, and Dt as a 
function of k for different restitution coefficients. 
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FIG. 2. Stationary temperature divided by the theoretical 
value as a function of the inelasticity 1 — a for different den- 
sities. The temperatures are computed in molecular dynamic 
simulations of the collisional model described by the collision 
rule ((35} while the theoretical values are obtained assuming 
Maxwellian distributions and absence of velocity correlations. 
From top to bottom p =0.1, 0.2, 0.3, 0.4, and p =0.8. Inset: 
Theoretical (continuous line) and simulational (symbols) di- 
mensionless stationary temperature T si /A 2 as a function of 
the restitution coefficient. The values for different densities 
collapse on the scale of the figure. 



A. Transport coefficients and thermodynamic 
properties 

The pseudo-hydrodynamic equations must be supple- 
mented by the equation of state p{p, T), energy injection 
rate G{p,T) and transport coefficients k, 77, and r\y . The 
computation of these require the analysis of the associ- 
ated Enskog equation to find first the stationary distribu- 
tion to compute p and G. Second, the Enskog equation 
must be analyzed using the Chapman-Enskog method to 
compute the transport coefficients. The purpose of this 
article is to analyze the hydrodynamic models of granular 
matter rather than performing a kinetic theory descrip- 
tion. There are numerous attempts to compute transport 
coefficients of granular fluids under different conditions, 
for example the homogeneous cooling state in two and 
three dimensions [Tol- [l2[ . the randomly driven gas (30j . 
the uniform shear flow ,[l5j. and others. It has become 
evident that the transport coefficients and equations of 
state depend strongly on the energy injection mechanism, 
the reference state and not only on the restitution coeffi- 
cient. Therefore, previous predictions of transport coeffi- 
cients or equations of state are not valid for the collisional 
model presented in this article. Even the first inelasticity 
correction is not valid. 

Considering the above discussion, in order to make 
quantitative comparison with the simulation results, we 
will use quasielastic values. That is, the transport coef- 
ficients and the equation of state are those of the elas- 
tic fluid and only the energy injection rate is computed 
considering the inelasticity Expressions for these 
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FIG. 3. Dynamic structure factor. The density and resti- 
tution coefficients are p = 0.4 and a = 0.94, respectively, 
while the wavevectors are ka = 0.01 (top) -corresponding to 
the dissipative regime- and ka = 0.06 (bottom) -quasielastic 
regime. 

functions are given in Appendix [B] 

Using the numerical value of the pseudo-hydrodynamic 
matrix at each wavevector, we perform a full diagonal- 
ization of it. By this procedure we obtain ujb, I\ and Z?t 
as a function of the wavevector, results that are com- 
pared with those obtained from the molecular dynamics 
simulations. 



B. Sound modes 

In the whole range of inelasticities the sound modes 
are visible, being possible to obtain lub and T with good 
accuracy. As predicted, the sound frequencies go linearly 
with k when the wavevectors are either in the inelastic 
or in the quasielastic regime, described in Sect. Ill and 
IV, respectively. Figure 2] shows the sound velocity ws/fc 
as a function of k for a series of restitution coefficients. 
At small k, the sound velocity takes a constant value 
that tends to the dissipative speed Cd while at large k 
it takes a different value that approaches the adiabatic 
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speed c s , recovering the inelastic and quasielastic regimes 
discussed in the previous section. 
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|B1 The agreement is excellent showing that both the hy- 
drodynamic model is appropriate and that the quasielas- 
tic transport coefficients give a good approximation of 
the dynamics of the model, at least in the range of in- 
elasticities presented in the figure. 

The crossover wavevector fco between the inelastic and 
quasielastic regimes is obtained fitting the sound velocity 
to a Lorentzian 



LO B /k 



Cd 



i + fc 2 /fc 2 ' 



(41) 



Figure [5] shows that fco is linear with the inelasticity I — a 
as predicted in Eq. ([32]l. A linear fit gives fc„ B = (0.530± 
0.005) (1 — a) j a. In the case of the sound damping, the 
crossover wavevector is obtained with a Lorentzian fit 
similar to (l4"Tj) except that in this case all coefficients are 
free. The result, shown in Fig. [5]is linear with 1 — a and 
a fit gives fcg = (0.525 ± 0.005) (1 — a) /a which coincides 
with the value obtained using ujb- Both results should be 
compared with the prediction (|34[) that, simplifies when 
G = to 



fco = 



r 2 n 4 



4c v p 2 pp p pT - 2p 2 p 2 r 



c 2 v P A P 2 p 



1/4 



G T . (42) 



Using the transport coefficients for quasielastic hard 
disks, we get a linear dependence with inelasticity fco = 
0.64(1 — a) I 'a, that compares well with the simulations. 
Note that the value does not need to agree as the proce- 
dures to obtain fco are not exactly equivalent. 



FIG. 4. Dimensionless sound velocity uis/k (top) and sound 
damping constant T/k 2 (bottom) as a function of the di- 
mensional wavenumber ka for different restitution coefficients 
(«=1.00 □, 0.99 •, 0.96 ■ 0.90A, 0.80 ♦ and 0.70 T). Points 
are the results of the simulations and the lines the theoretical 
predictions using the quasielastic transport coefficients. The 
solid horizontal lines are the predictions for an elastic fluid 
using Enskog theory: adiabatic velocity c s = 2.968 (top, solid 
line), dissipative (isothermal) velocity Cd = 2.044 (top, dotted 
line), and sound damping constant F = 2.954 (bottom). 

For small fc, the width V of the sound mode is propor- 
tional to fc 2 as predicted. Therefore, Fig. 0] shows T/k 2 
for a series of restitution coefficients. A similar crossover 
from the dissipative regime for small fc and the quasielas- 
tic regime for large fc is observed. At large values of fc the 
sound damping constant tends to the quasielastic value 
which is almost independent of the restitution coefficient, 
while for small fc, in the inelastic regime, the sound damp- 
ing constant diverges as 1/(1 — a). Note that the elastic 
limit is singular as the crossover wavevector vanishes as 
well. Therefore, as expected, no divergence is obtained 
in the elastic case. 

In Fig. [4] we present as solid lines the predictions of 
the hydrodynamic equations. They are obtained by full 
diagonalization of the hydrodynamic matrix M using the 
quasielastic transport coefficients described in Appendix 
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FIG. 5. Crossover wavevector for different inelasticities com- 
puted from Lorentzian fits of the sound speed (O) an d sound 
damping constant (•). The dotted line is the result of a linear 
fit, fc = (0.525 ± 0.005)(1 - a)/a. 



C. Heat mode 

In the quasielastic regime the heat mode is well defined 
and the fitting procedure gives accurate values for 7 and 
Dt- In the inelastic regime (high inelasticities and low 
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wavevectors), the heat mode can be hidden by the sound 
modes as shown in Fig. [3]and the dynamics is represented 
by two modes only. However, using the numerical data 
from the simulations it is possible to force a fit F(k,i) 
with three peaks using Eq. (fT9|) . Fig.[6]shows the relative 
amplitude of the heat peak compared to the sound peaks 
computed as the ratio between their areas 7 — 1. It is 
clear that in the inelastic regime the signal to noise ratio 
is poor and the precision in the fitted values for the heat 
mode is low. 



1.4 



0.20 



1.2 
1.0 
0.8 
^ 0.6 
0.4 
0.2 



n g n ° □ ° □ □ ° 



□ □ 



T-^ «-S-» 



• A f 



0.0 
0.00 



T ▼ 

T A A A A 



0.05 



0.10 

ka 



0.15 



0.20 



FIG. 6. Relative amplitude of the thermal peak to the sound 
peaks in the structure factor, Ir/Ib = 7 — 1, as a function 
of the dimensional wavenumber ka for different restitution 
coefficients (a=1.00 □, 0.98 •, 0.96 A, 0.92 ■ 0.88 v, 0.80 T 
and a =0.60 A). The solid horizontal line is the prediction 
for an elastic fluid using Enskog theory: 7 — 1 = 1.109. 

Figure [7] presents the fitted values of the width of the 
thermal peak, Dt- As expected, deep into the inelastic 
regime, it is not possible to obtain Dt with precision. In 
the quasielastic regime (large wavevectors) Dt shows a 
quadratic dependence with k. However, from this ten- 
dency it is not possible to obtain the homogeneous dissi- 
pation Gt by extrapolating it to k — > 0. This is due to an 
increase of Dt when decreasing k in the inelastic regime, 
as predicted in (|2"Tj) for which it can be verified that the 
coefficient of the k 2 term is negative (see Appendix |Bj . 
The position of the minimum of Dt is of the order of the 
crossover wavevector k^. 

Again, the comparison with the theoretical eigenvalues 
computed using the quasielastic transport coefficients is 
excellent for this range of inelasticities. 



D. Transverse mode 

The transverse dynamics is much simpler as it decou- 
ples from the longitudinal modes. As shown in Eq. (|9]) 
the transverse eigenvalue is Aj^ = k 2 v = k 2 r\j p. The 
transverse eigenvalue is obtained from the simulations 
computing the self-correlation function of the transverse 
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FIG. 7. Width of the termal peak Dt as a function of the di- 
mensional wavenumber ka for different restitution coefficients 
(a=1.00 □, 0.98 •, 0.96 ■ 0.92 ▲ and a = 0.88 T). Points 
are the results of the simulations and the lines the theoretical 
predictions using the quasielastic transport coefficients. 



current 



JV 



jx(k,t) = 2(i-kk).v j 



-ik-r ( (t) 



(43) 



where k = k/fc. The correlation function indeed decays 
exponentially, allowing the extraction of the transverse 
eigenvalue. In the range ka < 0.2 the transversal eigen- 
values are quadratic with k and the resulting viscosities 
are presented in Fig. [5] for different inelasticities. The 
elastic value agrees with the Enskog prediction in 2D and 
it takes smaller values as the inelasticity is increased. A 
linear fit gives 



v = (1.314 ± 0.004) - (0.37 ± 0.01)(1 - a) 



(44) 



in the range 0.1 < a < 1, that is almost up to the plastic 
limit. 



VIII. CONCLUSIONS 

We study the dynamics of a granular medium sub- 
jected to a bulk energy injection. Making simple 
generic assumptions on the injection method (momen- 
tum conservation and Galilean invariance) the pseudo- 
hydrodynamic equations are written. These describe the 
dynamics of the conserved density and velocity fields and 
the non-conserved temperature field. 

The fluctuations about the stationary state are ana- 
lyzed and described in terms of the dynamic structure 
factor and the corresponding eigenvalues of the inelastic 
hydrodynamic matrix. The dynamics near the stationary 
state is characterized in term of the following modes: the 
viscous mode that decouples as usual from the rest, the 
sound modes and the heat mode. Two regimes are clearly 
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FIG. 8. Dimensionless kinematic viscosity v = rj/(p/VT)&s a 
function of the inelasticity 1 — a obtained from simulations of 
the collisional model (symbols). The dashed line is a linear 
fit v= (1.314± 0.004) -(0.37±0.01)(1-q) and the solid line 
is the theoretical value for an elastic fluid obtained form the 
Enskog value Elastic = 1-303. 



distinguished. First the dissipative regime in which the 
heat mode is suppressed and the effective dynamics is 
reduced to only two fields and, second, the quasielastic 
regime in which the heat mode is visible. The crossover 
wavevector is proportional to the dissipation. In the dis- 
sipative regime the sound speed is isothermal and the 
sound damping constant becomes large, diverging in the 
limit of small dissipation and small wavevectors; the elas- 
tic limit is singular as it can be obtained by only making 
first the dissipation small and only later the wavevector 
can be small. In the quasielastic regime the sound veloc- 
ity is the adiabatic one and the sound damping constant 
is similar to the elastic value. 

The general predictions are compared with a collisional 
model that mimics the effective two-dimensional dynam- 
ics of a horizontally shallow three-dimensional system. 
In the shallow system the box is vertically vibrated and 
the energy is transferred from the box to the grains and 
later to the two-dimensional degrees of freedom through 
grain-grain collisions. To model this, we consider a purely 
two-dimensional granular fluid and, in collisions, a fixed 
additional separation velocity is added to the postcolli- 
sional velocities of each grain. A stationary temperature 
is reached that depends on the restitution coefficient and 
this added velocity. 

Molecular dynamics simulation of this model confirm 
the qualitative description of the hydrodynamic modes 
for the fluctuations. Besides, using the transport coeffi- 
cients of the elastic fluid plus the energy injection func- 
tion computed for the model, there is a very good nu- 
merical agreement with the theoretical predictions. It 
is difficult, however, to obtain the transport coefficients 
from a fit of the simulational eigenvalues. This is due to 
numerical accuracy and the fact that there are too many 
unknowns to be fit: the transport coefficients, the equa- 



tion of state and the energy injection rate. In the case 
of the transverse mode it is possible to fit the viscosity 
obtaining its dependence with the inelasticity 

ACKNOWLEDGMENTS 

The authors thank Ana Asenjo for useful com- 
ments. The research is partially supported by Spanish 
grants MODELICO and ENFASIS, the Fondecyt grants 
1100100 and 1120775, and Proyecto Anillo ACT 127. 

Appendix A: Dynamic Structure Factors 

In a system composed of N grains in a volume V 
(global density p = N/V), the local density field is de- 
fined as 

JV 

p(r, = (Al) 

i=l 

where i*i(t) is the position of the i-th particle at time t. 
The densiy-density correlation function is 

C pp (r, t) = (p(r + r',t + t')p(r', f )) - p 2 ■ (A2) 

For the velocity correlation function, which is a tensorial 
quantity, the dynamic variable is 

JV 

j(r ) t)=X)vi<y(r-r i (t)). (A3) 
i=l 

By taking the Fourier transform in space of Cab (r, t) 
we obtain the so called intermediate scattering function. 
For the density-density case, it reads 

F(k, t) = p- 1 J dr e- ik r C pp (r, t) (A4) 

= ^E ex PHk- (ri(t) -r,(0))]\ - (2nfpS(k). 

(A5) 

The last term containing a Dirac delta at k = only 
represents the mass conservation, that will be dropped 
from here on. 

Subsequently, one can perform a Fourier transform in 
time variable of F (fc, t) to obtain the dynamic structure 
factor, S{k,uS). Such structure factor is a fundamental 
tool in the study of fluid systems, like gases, liquids, 
polymers, and colloids [42[ . The reason is that, at long 
wavelengths, k — > 0, S(k,w) encodes many equilibrium 
and non equilibrium properties of the fluid. At long wave- 
lengths, it allows full evaluation, by the so called Landau- 
Plazceck approximation, where the evolution of the fields 
is given by the hydrodynamic equations. In this regime, 
S(k, ui) shows three Lorentzian peaks. Its expression is 
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26] (we quote it here for reference and comparison with 
the inelastic case) 



S(k,u)/S(k) 



1 



2D T k 2 



{D T k 2 ) 2 
Yk 2 



1 



-f(uj±c s k) 2 + (Tk 2 ) 2 ' 



(A6) 



The first term represents a peak located at cj = (called 
Rayleigh peak) and appears as a consequence of energy 
conservation. It has a width given by Dt = k / (pcy)k 2 , 
where n is the heat conductivity of the fluid, p is the 
average density and cy is the specific heat at constant 
volume. Such a peak carries the information about the 
entropy evolution in the system. 

There are two other peaks, Brillouin peaks (repre- 
sented by the symbol ± in the denominator), symmetric 
respect to that of uj = 0, located at uj — ±ic s k, where c s 
is the adiabatic speed of sound (c s = \J (c p /cv){dp/dp), 
where c p is the specific heat at constant pressure) . Their 
presence is a consequence of the conservation of momen- 
tum and the inertia of the system. The width of these 
peaks is Tk 2 , where T is the sound damping constant, 
a combination of shear and bulk viscosities and the heat 
conductivity. Moreover, the area enclosed by the thermal 
peak divided by the area under the sound peaks is 7 — 1, 
being 7 the adiabatic constant, so it yields the ratio of 
specific heats c p /c v . 

The expression at Eq. (|A6[) is obtained at the lowest 
order in the wave vector k. The next order term in k 
gives the so-called asymmetric contribution to the Bril- 
louin peaks that vanishes in the hydrodynamic limit, but 
have a finite contribution at finite k [261 ] . Such asym- 
metric peaks are considered in the paper. In summary, 
the measure of S(k,ui) in molecular fluids allows us to 
obtain transport coefficients and some thermodynamic 
properties. 



Appendix B: Thermodynamic properties and 
transport coefficients of the elastic hard disk fluid 

Here we provide the expressions for the thermody- 
namic properties and transport coefficients of the elastic 
hard disk fluid (2D), used in the comparison with the 
simulation results. Units are such that the grain diame- 
ter a and mass m are set to one. 

The thermodynamic properties of an elastic hard disk 
(2D) fluid can be obtained using the expressions of the 



equation of state [43, 44J and the specific heat at constant 
volume 



cy 



V = pT- 
1 



/ 2 2 

[ 1 4. EJL. 
\ 1 ' 128 



(Bl) 
(B2) 



where temperature is measured in energy units. From 
these it is straightforward to obtain the adiabatic and 
dissipative sound speeds, and the adiabatic constant 



7-1 



c s = ^Pp+P 2 T /cvP 2 « 2.97VT, 
p- p « 2.04VT, 



Cp/cy W 1.11, 



(B3) 

(B4) 
(B5) 



where, in the right equalities the reference density p = 0.4 
used in the simulations along the paper has been re- 
placed. 

The transport coefficients of the hard disk gas have 
been obtained using the Enskog transport equation (45j . 
and are valid up to moderate densities 



1.022 T 



V = 



T]V = 



X 

1.022 

X V 2vr 
1.029 Fxf 

X V 7T 



^ [1 + (2p X ) + 0.8729(2p X ) 2 ] 
[1.246(2 PX ) 2 ] a 0.26VT 



-(2 PX ) + 0.8718(2p X ) : 



0.52VT 
(B6) 

(B7) 

« 2.46^, 
(B8) 



where x is the pair correlation function at contact given 
by 



1 LZL£ 

x 16 4 



1.83, 



(B9) 



and, as before, the reference density p = 0.4 has been 
replaced in the right equalities. 

Finally, some expressions that are used in the compu- 
tations of the eigenvalues are 



D T = 



np p p 



p T + CvPpP 2 



2.92, 



2 



Kp 2 ^ 



2c v p{p T +p p p 2 c v ) 



2.95. 



(B10) 
(Bll) 
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